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Abstract 


Fighter  jets  and  other  aircraft  with  high  specific  thrust  engines  can  make  a  particularly  intense  type 
of  noise  that  has  come  to  be  called  crackle.  Its  rough,  sporadic  character  is  sometimes  likened  to 
super-loud  paper  tearing  or  the  sharp  staccato  of  water  drops  added  to  hot  oil.  It  is  unmistakably 
different  than  lower  speed  jet  noise.  Its  frequency  content,  extreme  intensity,  and  sporadic  character 
make  it  particularly  annoying  and  potentially  injurious  to  insufficiently  protected  hearing.  Observations 
suggest  that  it  consists  of  series  of  intense  shock-like  sound  waves,  which  arrive  intermittently  in  groups. 
Between  these  shock  groups  a  less  intense  but  still  loud  standard  jet  noise  is  heard.  The  steepened  shock¬ 
like  character  is  suggestive  of  nonlinear  development,  but  to  what  degree  they  attain  this  character  at 
generation  versus  nonlinear  steepening  during  propagation  is  unclear.  There  is  evidence  that  both  are 
key  factors  in  crackle.  It  is  also  unclear  whether  the  intermittency  of  crackle  is  due  to  the  changing 
character  of  the  turbulence  or  if  it  is  a  stochastic  propagation  effect. 

We  have  complete  a  set  of  detailed  direct  numerical  simulations  of  free  shear  flows  to  investigate 
turbulence  as  a  source  of  jet  crackle.  This  report  includes  (1)  detailed  documentation  of  these  simula¬ 
tions,  (2)  their  verification  and  validation,  (3)  and  several  post-processing  investigation  the  underlying 
mechanics  of  jet  crackle,  including  detailed  quantification  of  the  near  acoustic  field,  analysis  of  statistical 
metrics  of  jet  crackle,  analysis  of  sound-field  nonlinear  effects,  and  the  relation  of  crackle  to  the  linear 
instability  modes  supported  by  free  shear  flows.  This  is  the  first-ever  such  detailed  investigation  of  these 
mechanisms. 


1  Summary 

The  summary  of  our  accomplishments  of  this  funding  period  are: 

•  We  completed  a  detailed  study  of  the  turbulence  and  radiated  sound  of  high-speed  turbulent  shear 
layers  for  Mach  numbers,  M  =  0.9,  1.5,  2.5,  and  3.5.  The  results  indicate  that  for  M  >  2.5  crackle 
noise  correlated  to  presure  skewness  Skip')  >  0.4  is  present  in  the  very  near  acoustic  field.  This 
was  the  first  such  detailed  DNS  investigation  of  turbulence  as  a  source  of  crackle. 

•  We  confirmed  that  these  flows  have  broad-banded,  realistic  turbulence  and  that  low-order  statistics 
match  well  with  previous  simulation  and  their  planar  shear  layer  experimental  counterpart. 

•  Detailed  analysis  of  the  near-field  acoustics  shows  weak  shock-like  waves  undergo  nonlinear  inter¬ 
actions  during  sound  propagation  leading  to  elevated  skewness  in  the  very  near  field.  After  that, 
they  propagate  with  little  additional  nonlinear  distrotion. 

•  We  completed  a  higher-Reynolds  number  simulation  for  M  =  2.5  at  twice  the  grid  resolution. 
The  results  show  that  crackle  is  independent  of  Reynolds  number,  for  the  Reynolds  number  range 
considered. 

•  An  analysis  using  linear  stability  theory  of  a  parallel  shear  flow  was  completed.  It  is  clear  that  for 
M  >  1.2  (Mc  >  0.6)  the  flow  becomes  mor  e  three-dimensional  which  is  reflected  in  the  turbulence 
statistics. 

•  A  formulation  of  the  statistical  skewness  budget  was  evaluated  in  the  near  acoustic  field.  The 
results  indicate  statistical  skewness  originates  near  or  at  the  turbulent  source.  The  effect  of  heat 
conduction  and  viscosity  which  would  serve  to  decrease  statistical  skewness  is  nearly  balanced 
by  other  acoustic  contributions.  This  suggests  the  footprint  of  elevated  skewness  arises  near  the 
turbulence  and  propagates  nearly  unchanged  over  distances  considered  in  the  DNS. 

2  Introduction 

2.1  Background 

Military  aircraft  and  other  jet  engines  operating  with  high-specific  thrust  have  been  observed  to  radiate 
a  particularly  intense  and  distinct  sound  that  has  become  known  as  ‘crackle’.  Crackle  has  been  described 
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as  a  ‘startling  staccato  of  cracks  and  bangs’  (Ffowcs  Williams  et  al. ,  1975)  and,  although  a  very  annoying 
aspect  of  jet  noise  when  observed  from  a  distance,  this  distinct  sound  can  cause  aural  injury  for  personnel 
working  in  close  proximity.  Reducing  or  eliminating  crackle  from  the  overall  jet  noise  would  reduce 
annoyance,  improve  environmental  noise  levels  near  airbases,  and  decrease  the  potential  for  injury. 

Crackle  has  been  extensively  documented  after  it  was  first  reported  as  a  distinct  phenomenon  by 
Ffowcs  Williams  et  al.  (1975)  for  supersonic  full-scale  engines  and  lab-scale  jets.  Although  crackle 
is  an  extremely  intense  sound,  it  is  the  intermittency  of  crackle  that  makes  it  particularly  annoying 
(Ffowcs  Williams  et  al.,  1975).  At  seemingly  random  times,  crackle  lasts  between  0.1  and  1  seconds  at  1 
to  2  second  intervals  (Ffowcs  Williams  et  al.,  1975).  During  this  period,  strong  shock-like  compressions 
are  observed  at  a  frequency  of  10~2  seconds.  The  time  trace  of  acoustic  pressure  during  a  crackle  event 
reveals  a  peculiar  feature — strong,  steepened  compressions  (lasting  1  millisecond)  followed  by  longer 
weaker  expansions  (lasting  2  to  3  milliseconds)  with  the  expansion  amplitude  about  |  the  compression 
peak.  The  authors  noticed  that  the  statistical  skewness  of  acoustic  pressure  would  quantify  the  amplitude 
asymmetry  between  the  sharp,  intense  peaks  and  rounded,  weaker  valleys  in  the  sound  signal.  Correlating 
the  perception  of  crackle  with  statistical  skewness, 

Sk(p)  2^3/2’  W 

It  was  observed  that  sound  with  Skip1)  >  0.4  would  be  noticeably  crackling.  Though  not  a  rigorous  defini¬ 
tion,  statistical  skewness  is  a  suitable  metric  for  correlating  to  the  perception  of  crackle  (Ffowcs  Williams 
et  al,  1975;  Krothapalli  et  al. ,  2000;  Petitjean  &  McLaughlin,  2003;  Petitjean  et  al.,  2005).  Gee  et  al. 
(2007)  suggests  that  the  skewness  of  the  time  derivative  of  pressure  is  a  suitable  metric  for  quantifying 
sharp  changes  in  the  pressure  signal,  which  are  also  known  to  occur  in  crackle  noise.  Other  metrics  such 
as  pressure  spectra  are  insensitive  to  detecting  the  presence  of  crackle  (Ffowcs  Williams  et  al.,  1975)  in 
part  because  of  its  intermittency:  crackle  events  occurring  in  about  5%  of  the  acoustic  signal  (Krothapalli 
et  al,  2000). 

Crackle  levels  depend  on  jet  operating  conditions.  Skip1)  is  known  to  increase  with  increasing  flow 
velocity  (Ffowcs  Williams  et  al,  1975;  Szewczyk,  1978;  Krothapalli  et  al.,  2000;  Petitjean  &  McLaughlin, 
2003;  Petitjean  et  al,  2005;  Gee  et  al,  2013)  and  crackle  is  apparent  when  the  jet  velocity  is  high 
enough  that  turbulent  eddies  in  the  shear  layer  move  supersonically  relative  to  the  ambient  free  stream 
(Krothapalli  et  al,  2003).  Initially  thought  to  be  insensitive  to  jet  temperature  (Ffowcs  Williams  et  al, 
1975;  Szewczyk,  1978),  within  the  temperature  range  tested,  crackle  levels  do  correlate  with  higher 
temperature  ratios  (Krothapalli  et  al,  2000;  Petitjean  &  McLaughlin,  2003).  The  frequency  of  the 
sporadic  crackle  events  also  increases  with  increasing  jet  velocity  and  temperature.  Crackle  does  not 
seem  to  be  correlated  to  the  presence  of  internal  shock-cells  from  imperfect  expansion  (Ffowcs  Williams 
et  al,  1975;  Szewczyk,  1978;  Krothapalli  et  al. ,  2000)  nor  to  rough  combustion  processes  present  in  full- 
scale  engine  tests  (Ffowcs  Williams  et  al,  1975).  High  skewness  levels  have  been  correlated  to  angles  near 
the  expected  Mach  angle  (Ffowcs  Williams  et  al,  1975;  Laufer  et  al,  1976;  Szewczyk,  1978;  Krothapalli 
et  al,  2000;  Petitjean  &  McLaughlin,  2003;  Petitjean  et  al,  2005;  Gee  et  al,  2013)  suggesting  that  crackle 
may  be  associated  with  the  Mach  wave  radiation  mechanism  known  to  occur  in  supersonic  free-shear-flow 
turbulence,  though  it  has  been  observed  that  flows  radiating  Mach  waves  may  not  radiate  noise  perceived 
to  crackle  (Krothapalli  et  al,  2000). 

It  appears  that  crackle  is  correlated  to  signals  with  intermittent  periods  of  steepened  shock-like  waves 
followed  by  weaker,  longer,  rounded  rarefaction  regions,  but  to  what  extent  this  character  originates  at 
the  turbulent  source  or  arises  from  nonlinear  propagation  effects  is  unclear.  Ffowcs  Williams  et  al  (1975) 
suggested  that  crackle  originates  at  or  near  the  source.  Within  the  range  that  the  acoustic  measurements 
were  made,  Ffowcs  Williams  et  al  (1975)  demonstrated  that  crackle  was  independent  of  jet  scale  and 
propagation  distance,  which  suggests  that  cumulative  nonlinear  propagation  effects  did  not  distort  the 
signal.  Visualizations  do  suggest  that  the  radiated  waves  are  steepened  near  the  shear  layer  of  the 
jet  (Krothapalli  et  al.,  2000;  Freund  et  al.,  2000;  Darke  &  Freund,  2001;  Nichols  et  al,  2013)  but  the 
connection  between  the  steepened  signal  footprint  near  the  source  and  the  signal  observed  in  the  acoustic 
field  has  not  been  addressed  in  detail.  There  is  evidence  that  Sk(p')  does  evolve  with  propagation  distance 
(Szewczyk,  1978;  Petitjean  &  McLaughlin,  2003;  Petitjean  et  al,  2005)  with  noticeable  nonlinearity 
during  its  course.  The  variation  of  sound  spectra  with  propagation  distance  indicates  energy  being 
transferred  from  mid-frequencies  to  low-  and  high-frequencies  (Petitjean  &  McLaughlin,  2003;  Petitjean 
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et  al,  2005;  Fievet  et  al.,  2014).  The  expectation  is  that  wave  merging  due  to  nonlinearity  will  increase 
low-frequency  energy  content,  though  there  is  no  direct  observation  of  this,  and  wave-steepening  will 
contribute  energy  to  high-frequency  modes. 

Weak-shock  theory,  when  applied  to  the  propagation  of  broadband  noise,  predicts  a  final  wavetrain  of 
N-waves  (Lighthill,  1956;  Pestorius  &  Blackstock,  1974;  Freund  et  al.,  2000)  with  symmetric  amplitude 
distribution  (i.e.  Skip')  =  0),  not  the  asymmetric,  peaked  signals  linked  to  the  perception  of  crackle  noise. 
This  suggests  that  skewness  arises  from  some  other  factor  than  one-dimensional  weak-shock  propagation. 
Crighton  (1986)  showed  that  linear  dispersion  and  attenuation  coupled  with  nonlinear  mechanisms  will 
increase  Sk  with  distance,  which  might  explain  the  skewness  variations  observed  by  Szewczyk  (1978), 
Petitjean  &  McLaughlin  (2003),  and  Petitjean  et  al.  (2005). 

2.2  Objectives  and  outline 

The  present  simulations  are  designed  to  provide  a  detailed  description  of  the  turbulence  mechanisms 
and  any  near-Held  nonlinearity  that  leads  to  acoustic  fields  of  the  kind  associated  with  the  perception 
of  crackle.  The  simulations  resolve  all  the  energetic  scales  with  no  sub-grid  scale  models,  which  would 
potentially  affect  important  features  of  the  flow  field,  such  as  the  sharp  shock-like  compressions  in  the 
near  acoustic  field.  Though  the  shock  thickness  is  set  by  the  viscosity  for  the  parameters  considered, 
there  is  a  small,  but  finite  numerical  dissipation  in  the  simulations  which  is  expected  to  weakly  thicken 
shocks.  While  jet  noise  is  the  most  important  engineering  situation,  we  chose  to  study  a  temporally 
developing  planar  shear  layer  because  it  provides  a  clearer  perspective  of  the  root  mechanisms  of  this 
sound  generation. 

Though  the  equations  are  frame  invariant,  some  care  is  required  in  switching  intuitively  from  the 
supersonic  advection  generating  Mach-like  waves  in  a  frame  attached  to  a  nozzle  to  the  corresponding 
phenomenon  in  a  time-developing  flow,  where  transient  approximately  stationary  (on  average)  turbulence 
features  are  adjacent  to  a  supersonic  stream.  This  configuration  can  be  considered  a  model  for  the  near¬ 
nozzle  region  of  a  high-Reynolds-number  jet,  where  the  turbulence  is  concentrated  in  a  weakly  curved  (for 
a  typical  round  jet)  shear  layer  between  the  high-speed  potential  core  flow  and  the  surrounding  co-flow. 
Because  the  simulations  focus  on  this  small  region,  it  does  allow  us  to  explicitly  resolve  a  broader  range 
of  turbulence  scales  than  could  be  represented  in  a  full  jet  simulation.  In  essence,  it  provides  a  Reynolds 
number  realistic  representation  of  a  piece  of  a  high-Reynolds-number  jet,  and  it  therefore  allows  us  to 
probe  the  sound  generation  mechanisms  of  high-Reynolds-number  turbulence.  It  can  be  anticipated  the 
acoustic  far-field  will  not  include  the  geometric  propagation  effect  associated  with  a  round  jet,  though 
this  will  not  be  a  large  effect  for  studying  the  very  near  acoustic  field.  The  flow  configuration  is  not 
ideally  suited  for  making  point-to-point  comparisons  with  lab  measurements  because  acoustic  data  is 
typically  acquired  at  >  10 Dj  from  the  shear  layer,  much  farther  away  than  our  simulations,  and  there 
is  no  direct  experimental  counterpart  to  temporally  developing  shear  layers.  We  do  expect  that  lab 
measurements  would  ‘see’  a  steepened,  skewed  footprint  of  a  weakly  nonlinear  signal.  The  nearly  linear 
signal  would  continue  to  propagate  with  its  peculiar,  asymmetric  character  and  undergo  slow,  cumulative 
nonlinear  changes  associated  with  propagation  of  this  type.  Near-field  acoustic  measurements  support 
this  hypothesis  (Petitjean  &  McLaughlin,  2003;  Petitjean  et  al.,  2005;  Baars  &  Tinney,  2014;  Fievet 
et  al.,  2014).  The  details  related  to  the  simulation  configurations  and  numerical  methods  are  provided 
in  section  3. 

Very  near-field  sound  of  the  temporally  developing  flow  does  indeed  exhibit  a  sound  field  with  the 
hallmarks  of  crackle:  Sk{p')  >  0.4  with  steepened,  compressions  in  the  sound  field.  Analysis  of  our 
database  will  show  the  key  features  of  the  near-field  acoustics  that  lead  to  these  characteristics.  Key 
results  of  our  simulations  including  shear  layer  characteristics,  turbulence  statistics,  and  the  evolution  of 
the  Mach-like  wavefielcl  leading  to  the  sound  correlated  to  crackle  is  given  in  section  4.  A  supplementary 
lrigher-Reynolds  number  simulation  for  M  =  2.5  shows  that  crackle  is  insensitive  to  Reynolds  number. 

In  section  5  a  formulation  of  skewness  budget  is  given.  The  budgets  for  a  nominally  non-crackling 
and  crackling,  M  =  0.9  and  2.5,  respectively,  are  compared  and  we  show  that  for  the  low-speed  case,  the 
entropic  contribution  to  skewness  is  small  (no  shocks  are  present)  and  that  skewness  remains  unchanged 
because  acoustic  mechanisms  integrate  roughly  to  zero.  On  the  other  hand,  for  M  =  2.5,  weak  shocks 
are  present  in  the  acoustic  field  increasing  the  entropic  sink  of  skewness.  This  sink  is  balanced  by  other 
acoustic  mechanisms.  The  skewness  here  also  remains  unchanged  suggesting  the  footprint  of  skewness, 
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a  metric  correlated  to  crackle,  occurs  at  or  very  near  the  turbulent  source. 


3  Simulation  details 

3.1  Configuration 

The  three-dimensional  compressible  flow  equations  in  Cartesian  coordinates  were  solved  for  the  tempo¬ 
rally  developing  shear  layer  shown  in  figure  1.  The  details  of  the  full  governing  equations  are  provided 
in  appendix  A.  The  flow  domain  dimensions  are  LxxLyxLz  =  120  8^  x  200  80  <5^,  where  S ^  is  the 

initial  momentum  thickness, 


'•Ly/2 


5m{t)  = 


PooAU2 


p(\AU  —  u)(hAU  +  u)  dy , 


(2) 


'-Ly/2 

at  time  t  =  0.  In  (2),  ()  and  ()  indicate  Reynolds  and  Favre  averages,  respectively.  The  Reynolds  number 

pocAUSrn(t') 


= 


(3) 


was  Rtsm  =  60  initially  and  reached  2100  when  =  35.  A  supplementary  simulation  at  an  initial 

Reynolds  number  of  Res^  =  120  for  M  =  2.5  was  simulated  until  Regm  =  4000.  Details  regarding  this 
simulation  are  given  in  appendix  C. 

The  planar  shear  layer  was  initialized  with  divergence-free  velocity  perturbations  using  the  approach 
of  Kleinman  &  Freund  (2008)  added  to  a  streamwise  velocity  profile  given  by 


A  U  (  v  \ 
u(y)  =  ~y~  tanh  J  ' 

Four  flow  speeds — M  =  0.9,  1.5,  2.5,  and  3.5 — were  simulated  until  the  momentum  thickness  reached 
35  times  the  initial  momentum  thickness.  The  Mach  number  is  defined  as  M  =  AU/Coo ,  where  AC  is 
the  velocity  difference  between  the  free-streams.  The  M  =  0.9,  1.5,  and  3.5  cases  were  initialized  by 
rescaling  the  M  =  2.5  velocity  field  to  their  respective  freestream  velocity  values. 


Figure  1:  Simulation  schematic. 
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3.2  Numerical  methods 

This  domain  was  discretized  with  Nx  x  Ny  x  Nz  =  1536  x  2048  x  512  mesh  points.  The  solution  was 
integrated  in  time  using  a  standard  fourth-order  Runge-Kutta  method.  Derivatives  in  the  periodic 
streamwise  (x)  and  cross-stream  (y)  directions  are  computed  with  resolution-optimized  finite  differences. 
For  the  ^-derivative,  the  fourth-order,  dispersion-relation-preserving  scheme  of  Tam  &  Webb  (1993)  was 
used.  Though  Fourier  methods  could  have  been  used  in  the  x-periodic  direction,  a  local  finite-difference 
method  was  chosen  to  be  more  computationally  efficient  than  the  global  Fourier  method  since  the  domain 
was  parallelized  in  this  direction.  The  y-derivative  was  computed  with  an  eighth-order  pentadiagonal 
compact  finite-difference  scheme  (Lele,  1992).  Both  x-  and  y-second  derivatives  were  approximated  using 
a  maximum-order  explicit  sixth-order  scheme.  Near-boundary  grid  points  in  the  y-direction  had  reduced- 
order  central  differences  and  one-sided  differences  at  the  domain  extent.  It  is  assumed  that  the  low-order 
boundary  closure  methods  for  the  derivatives  are  not  expected  to  degrade  the  solution  since  they  are  far 
from  the  turbulence  and  inside  absorbing  sponge  zones  whose  effect  will  be  discussed  later.  The  periodic 
spanwise  ^-direction  used  Fourier  methods  for  both  first  and  second  derivatives.  Details  with  regard  to 
each  numerical  method  are  found  in  appendix  B. 

Fourier  analysis  of  the  finite-difference  approximations  provides  an  approximate  dispersion  relation 
to  compare  to  the  exact  relation.  Figure  2  and  3  shows  the  resolution  properties  of  the  first  and  second 
derivatives,  respectively. 


- exact; - cc-dir.;  -  -  -y-dir.;  . z-dir. 


Figure  2:  Modified  wave  number  (fcA)'  resolution  of  first-order  derivatives. 


Because  of  finite  computational  resources,  the  domain  in  the  y-direction  must  be  truncated,  and  it 
is  desired  that  the  computational  boundary  mimics  the  properties  of  a  radiation  boundary  conditions 
(i.e.  no  numerical  reflections  into  the  domain  having  a  harmful  effect  on  the  solution).  Experience 
has  shown  for  flows  of  this  kind  (Kleinman  &  Freund,  2008)  that  a  combination  of  an  aborbing  sponge 
region  (Freund,  1997&)  with  one-dimensional  characteristic  radiation  boundary  conditions  (Thompson, 
1987)  can  achieve  the  radiation-like  behavior  provided  the  amplitude  of  fluctuating  quantities  leaving 
the  domain  is  sufficienctly  small.  Sponge  regions  of  width  ws  =  44  <5^  are  modeled  at  y  =  ±Ly/2  by 
adding  a  source  term  to  the  right-hand-side  of  the  governing  equations  of  the  form 


C(2/)(Q  Qtarget)- 

The  target  solution  vector  is  defined  as 


Qtarget  — 


Poo  i  Poo 


Poo 

7“  1 


+  2  Poo 


(5) 

(6) 
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- exact;  -  -  -  x-  &  y-dir.; 


z-dir. 


kA 


Figure  3:  Modified  wave  number  (kA)"  resolution  of  second-order  derivatives. 


which  forces  the  solution  vector  q  to  the  desired  ambient  boundary  conditions.  The  sponge  strength  is 
quadratically  increasing  toward  the  domain  extent  given  by 


f. 

1 

4? 

1 

2 

i  \y\  >  if  -  ws 

II 

2 

ws 

(7) 

[o, 

otherwise 

To  suppress  the  mild  instabilities  these  numerical  schemes  are  known  to  have,  high-resolution  nu¬ 
merical  filtering  was  applied  in  all  three  coordinate  directions,  as  is  often  done  (Kleinman  &  Freund, 
2008).  In  the  spanwise  direction,  the  top  20%  of  wave  numbers  of  the  Fourier  transformed  solution 
were  filtered  to  zero.  Explicit  and  implicit  filtering  as  designed  by  Lele  (1992)  was  used  in  the  x-  and 
y-directions,  respectively.  In  x,  the  coefficients  of  the  fourth-order  scheme  were  chosen  so  the  transfer 
function,  equation  (C.2.2)  in  Lele  (1992),  is  T(k)  =  0.99  at  wavenumber  k  =  18^.  In  y1  the  coefficients 
of  the  sixth-order  scheme  resulted  T(k)  =  0.95  at  k  =  25^  and  T(k)  =  0.5  at  k  =  2.6  5^.  The  filtering 
was  done  on  the  entire  solution  at  each  time  step.  Numerical  tests  which  will  explained  later  show  that 
the  filtering  affects  the  high-wave-number  range  of  the  flow  which  contains  little  energy.  The  largest 
scales  are  unaffected  by  this  stabilizing  technique. 


- exact; - rc-dii. ;  -y-dir.;  . z-dir. 


<1 

-ie 


kA 


Figure  4:  Transfer  function  T(kA)  of  selective  high-wavenumber  filtering  schemes. 
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4  Results 

4.1  Turbulence  statistics 

The  growth  of  turbulent  shear  layers  shown  in  figure  5.  Though  there  is  an  initial  transient,  which  is 
known  to  occur  for  simulations  of  this  kind  (Pantano  &  Sarkar,  2002;  Kleinman  &  Freund,  2008),  the 
growth  appears  approximately  linear  compared  to  the  regression  between  20  <  <  35. 


- M  =  0.9;  ---  1.5;  - 2.5;  3.5 


Figure  5:  Shear  layer  growth  with  linear  regression  for  >  20. 


Taking  the  time  derivative  of  equation  2  and  neglecting  molecular  dissipation  effects  on  the  mean  flow, 
the  momentum  thickness  growth  rate  is 


Sm(t) 


f*Ly  /2 


(8) 


where  (”)  represents  the  fluctuation  by  subtracting  the  Favre  average  from  the  total.  Figure  6  shows  the 
growth  rate  of  the  simulations  compares  well  with  spatially  developing  experiments  (Elliott  &  Samirny, 
1990;  Goebel  &  Dutton,  1991;  Debisschop  et  al. ,  1994)  and  simulation  (Pantano  &  Sarkar,  2002;  Kleinman 
&  Freund,  2008).  As  expected,  shear  layer  growth  is  suppressed  with  increasing  M  and  compares  well 
with  previous  results. 


^  present  DNS 
□  Kleinman  &  Freund  (2008) 

•  Pantano  &  Sarkar  (2002) 
x  Elliott  &  Samirny  (1990) 

*  Goebel  &  Dutton  (1991) 
o  Debisschop  et  al.  (1994) 

0  0.5  1  1.5  2 

Mc 

Figure  6:  Shear  layer  growth  rate. 
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Table  1  compares  spatial  resolution  of  previous  compressible  shear  layers  simulations.  For  the  present 
simluations,  the  resolution  coupled  with  a  sufficiently  low  Reynolds  number  is  expected  to  be  adequately 
resolved. 


Reference 

Mc 

ReS"n 

A  x/5°m 

Aj/min /Sm 

A  z/50m 

Pantano  &  Sarkar  (2002) 

0.3-0. 7 

160 

0.672 

0.672 

0.672 

Pantano  &  Sarkar  (2002) 

1.1 

160 

0.675 

0.675 

0.677 

Kleinman  &  Freund  (2008) 

0.9 

69 

2.94 

2.34 

4.49 

Kleinman  &  Freund  (2008) 

0.9 

207 

0.976 

0.79 

1.47 

Kleinman  &  Freund  (2008) 

0.9 

414 

0.976 

0.79 

1.47 

Zhou  et  al.  (2012) 

0.7 

100 

0.673 

0.775 

0.631 

present  DNS 

0.9-3. 5 

60 

0.94 

1.09 

1.88 

present  DNS 

2.5 

120 

0.5 

0.5 

0.5 

Table  1:  Comparison  of  resolution. 

Figure  7  shows  the  spectrum  of  velocity  fluctuations  at  y  =  0,  the  midplane  of  the  mixing  layer.  The 
turbulence  is  broadbanded  and  appears  realistic. 
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Figure  7:  One-dimensional  spectra  of  it- velocity  fluctuations  at  y  =  0  when  =  35  in  the  (a) 

streamwise  and  (b)  spanwise  direction. 


Figure  8  shows  the  effect  the  numerical  stabilization  method  has  on  the  highest  wavenumbers.  There 
is  no  apparent  effect  on  the  largest,  energy-containing  scales. 


The  peak  turbulence  intensities  are  compared  with  previous  experiment  (Elliott  &  Samimy,  f990; 
Goebel  &  Dutton,  1991;  Debisschop  et  al. ,  1994)  and  simulation  (Pantano  &  Sarkar,  2002;  Kleinman  & 
Freund,  2008)  in  figure  9.  As  expected,  the  basic  trend  of  the  intensities  decreases  with  increasing  Me. 
Average  Reynolds  stress  profiles  esemble  averaged  in  time  during  approximate  linear  growth  between 
<5m(f)/<5m  =  20  and  35  are  shown  in  figure  10. 
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Figure  8:  One-dimensional  spectra  of  u-velocity  fluctuations  at  y  =  0  when  (5m(f)/<5^  =  35  for  M  =  2.5 
in  the  (a)  streamwise  and  (b)  spanwise  direction.  The  comparison  between  the  number  of  time  steps 
between  filter  applications  is  made  after  one  flow-through  time,  tfi  =  1440  S^/AU. 


The  Reynolds  stress  profiles  for  M  =  0.9-3. 5  are  shown  in  figures  11-14.  For  higher  Mach  numbers 
M  >  1.5,  the  self-similar  collapse  of  the  profiles  is  not  as  obvious  for  the  M  =  0.9. 


4.2  Near-field  acoustics 

Whether  or  not  Skip1)  is  a  useful  metric  for  the  perception  of  crackle,  it  is  undoubtedly  a  peculiar 
feature  of  the  sound  held.  We  start  by  confirming  that  the  simulations  do  produce  this  hallmark  feature 
of  nominally  crackling  jets.  Figure  15  shows  that  Skip')  at  all  y  increases  with  M.  This  appears  to  happen 
continuously,  without  a  distinct  transition  to  finite  Sk  and  nominally  ‘crackling’  behavior,  though  the 
discrete  M  simulated  does  not  preclude  this.  For  M  >  2.5,  the  Skip')  exceeds  0.4  expected  for  turbulent 
jets.  All  statistics  presented  throughout  this  section  have  been  averaged  over  the  x-  and  ^-directions  and 
ensemble  averaged  in  time  between  tAU/Sm{t)  =  1500  and  2500. 

For  any  M,  into  the  very  near  acoustic  held,  we  observe  an  approximately  constant  Skip')-  This 
may  seem  to  contradict  the  evidence  that  skewness  increases  with  propagation  distance  (Szewczyk,  1978; 
Petitjean  &  McLaughlin,  2003;  Petitjean  et  al,  2005),  though  signihcant  variation  is  not  expected  for 
the  propagation  distances  we  consider  here. 


Figure  16  shows  typical  pressure  traces  at  different  Mach  numbers.  These  have  sharper  compressive 
peaks  than  expansions  for  the  highest  speeds,  especially  the  M  =  3.5  case,  which  is  comparable  to 
experimental  observation  (Ffowcs  Williams  et  al.,  1975;  Laufer  et  al.,  1976;  Szewczyk,  1978;  Krothapalli 
et  al.,  2000;  Baars  &  Tinney,  2014).  Note  that  these  pressure  traces  are  close  to  the  turbulence,  at 
y  =  5.5  5m.  Thus,  skewed  pressure  signals  and  apparently  steepened  wave  prohles  exist  very  near  the 
turbulence. 
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Figure  9:  Peak  turbulent  intensities. 


4.3  Near-field  nonlinearity 

Directional  waves  radiating  from  high-speed  free-shear  flows  have  been  observed  by  schlieren  and  shad¬ 
owgraph  techniques  for  many  years.  The  orientation  of  these  directional  waves  provide  a  clue  to  the 
relative  speed  of  the  turbulence  source  that  generates  this  Mach-like  wave.  Laufer  (1961)  used  hot-wire 
anemometer  techniques  to  measure  u'-velocity  perturbations  in  the  acoustic  field.  His  measurements  also 
indicated  that  the  fluctuations  away  from  the  turbulent  boundary  layer  were  soley  acoustic  ones.  Duan 
et  al.  (2014)  also  provided  correlations  in  the  near  acoustic  field  that  suggest  the  Mach- like  waves  are 
backward- facing  waves  with  ( u'p ')  «  —  1  and  that  they  are  mainly  acoustic,  with  negligible  entropic  and 
vortical  contributions.  Assuming  a  purely  acoustic  field  with  planar  waves,  Laufer  (1961)  used  the  linear 
relation  between  pressure  and  the  particle  velocity  u'n  given  here  by 


<  =  J_P^ 

u  M 

and  with  the  geometric  relation  between  the  horziontal  and  normal  velocity  perturbation 


(9) 


v!  =  u'n  cos(6),  (10) 

the  orientation  of  wave  can  be  calcualted.  The  inferred  angle  from  the  measurements  of  u'  agreed  visually 
with  the  observed  wave  orientation  in  the  near  acoustic  field(Laufer,  1961).  Simulation  of  radiated  sound 
from  high-speed  turbulent  boundary  layers  shows  a  similar  qualitative  agreement  using  the  same  analysis 
(Duan  et  al .,  2014). 

Using  a  time-delay  cross-correlation  of  the  pressure  field  between  two  measurement  points,  Nichols 
et  al.  (2013)  indirectly  measured  the  angle  of  propagation,  and  the  orientation  qualitatively  agrees  with 
the  orientation  of  the  wave-like  structures  in  the  sound  field. 
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Figure  10:  Average  Reynolds  stress  profiles. 


Our  approach  measures  the  angle  of  structures  directly  from  the  flow  solution.  Figure  17  shows  the 
method  for  computing  wave  angle  at  a  single  point  in  the  near  acoustic  held.  The  angle  is  given  by 


9  =  cos  1 n 


yi 


(11) 


where  ny  is  the  y-component  of  the  normal  vector  computed  from  the  gradient  of  the  pressure  held 

Vp 


|Vp| ' 


(12) 


Wave  angle  statistics  are  chosen  to  be  conditional  on  the  strength  of  the  compression  specifically  a 
dilatation  threshold.  Figures  18-20  show  the  average  wave  angle  conditioned  on  dilatation  threshold  for 
M  =  1.5,  2.5  and  3.5.  The  M  =  0.9  case  was  omitted  from  these  results  because  Mach  wave  radiation  is 
not  expected  to  occur  and  the  sound  held  for  M  =  0.9  did  not  have  sufficient  compressions  for  the  range 
of  dilatation  thresholds  considered  to  accrue  converged  statistics.  Based  on  the  proposed  mechanism  of 
Mach  wave  radiation  from  high-speed  turbulent  shear  layers  (Phillips,  1960;  Ffowcs  Williams  &  Maidanik, 
1965),  it  is  expected  that  the  mean  wave  angle  in  the  near  acoustic  held  decreases  with  increasing  how 
speed:  (9)  =  55°,  34°,  and  27°  at  y  =  15<5m(f)  for  M  =  1.5,  2.5,  and  3.5,  respectively. 

The  range  of  wave  angles  are  quantihed  by  the  standard  deviations  in  figures  18-20.  The  width  of  one 
standard  deviation  suggests  a  range  of  source  speeds  within  the  turbulence.  Smaller  standard  deviations 
about  the  mean  are  observed  for  slower  speeds.  For  M  =  3.5,  in  particular,  the  range  of  angles  decreases 
with  distance  from  the  mixing  layer.  This  would  be  consistent  with  the  geometric  battening  and  possible 
merging  of  different  waves  due  to  weak-shock  propagation  nonlinearity. 

Figure  21  shows  the  effect  of  varying  the  dilatation  threshold  on  the  average  wave  angle.  For  M  >  2.5, 
decreasing  dc  which  includes  stronger  waves  suggests  the  trend  of  decreasing  wave  angle  with  distance 
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Figure  11:  Reynolds  stress  profiles  for  M  =  0.9  at  different  shear  layer  thicknesses. 


is  converging.  For  M  =  1.5,  the  trend  with  distance  seems  to  diverge  primarily  because  less  waves  are 
being  counted  which  is  later  confirmed  in  figure  23. 


Cumulative  propagation  effects  are  anticipated  to  reduce  the  number  of  distinct  waves  (Lighthill, 
1956).  Described  as  ‘wave-bunching’,  stronger  shocks  with  a  higher  propagation  velocity  will  ‘consume’ 
weaker,  more  slowly  propagating  ones  (Lighthill,  1956,  1993).  To  quantify  this  wave-merging  effect,  we 
sum  the  number  (IV)  of  contiguous  compression  regions  in  the  ^-direction  that  satisfies  the  threshold 
condition,  V  •  u  <  dc.  We  see  in  figure  22  a  qualitatively  different  evolution  of  wave  density  between 
the  M  =  1.5  and  the  nominally  crackling  cases.  The  wave  density  in  y  increases  for  M  =  1.5  because 
at  earlier  times  the  smaller  scales  of  the  turbulence  produce  more  distinct  waves.  At  larger  y  these 
shorter  waves  (and  higher  frequencies)  persist.  Nonlinear  merging  counters  this.  For  M  =  2.5  and 
3.5,  the  number  of  waves  decreases  into  the  acoustic  field.  The  effect  of  varying  dc  between  —0.1  and 
— 0.003<5m(f)/A/7  causes  the  NSm/Lx  curves  to  shift  downward  as  fewer  waves  are  being  counted,  though 
the  shape  remains  unchanged. 

The  effect  of  dilatation  treshold  is  again  considered  for  the  wave  counting  as  shown  in  figure  23.  Over¬ 
all,  the  basic  trends  with  distance  from  the  centerline  is  captured  regardless  of  the  dilatation  threshold. 
The  number  of  waves  counted  shifts  vertically  with  decreasing  dc  because  fewer  waves  are  being  averaged. 
For  the  smallest  dilatation  threshold  in  M  =  0.9,  the  trend  with  distance  is  not  preserved  because  the 
number  of  waves  contributing  to  the  mean  is  so  small. 
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1.5  at  different  shear  layer  thicknesses. 


The  statistics  of  mean  wave  angle  and  scaled  wave  density  suggest  waveforms  undergo  nonlinear 
interactions  in  the  near  acoustic  field.  Visualization  and  tracking  of  individual  wave  fronts  confirm  that 
waveforms  nonlinearly  interact  and  merge  into  a  single  wave  front  in  figure  24. 


The  two-point  correlation  of  fluctuating  pressure, 

Cp'p'  (S,  t)  =  (p'(x  +  <My(x,i)},  (13) 

provides  a  measure  of  the  mean  wave  structure.  We  consider  a  two-dimensional  correlation  in  which 
5  =  (5X,  5y .  0).  The  integral  length  scale, 

y  J  Cp>p>(6x,6xtan(£),0,t)d5x,  (14) 

along  a  line  inclined  at  angle  £  is  shown  in  figure  26.  For  M  =  3.5  case  the  length  of  the  wave  grows 
with  increasing  distance  from  the  turbulence  (figure  26)  as  would  be  consistent  for  merging  waves. 
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5  Pressure  variance  and  skewness  budget 

5.1  Derivation 

Thermodynamic  quantities  are  generally  expressed  as  a  function  of  any  two  thermodynamic  quantities. 
The  pressure  is  written  as  a  function  of  density  and  entropy 


p  =  p{p,s), 

and  it  follows  that 


(15) 


Dp  / dp\  Dp  / dp\  DS 
Dt  \  dp  )s  Dt  \  ds  )p  Dt  ’ 

where  the  operator  D / Dt  is  defined  as 

D  _  d  d 
Dt  dt  Ul  dxi 

Assuming  an  ideal  gas,  the  thermodynamic  relations  are  given  as 


(16) 


(17) 


(18) 


where  cv  is  the  specific  heat  at  constant  volume  and  c  is  the  isentropic  speed  of  sound  given  as 
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Figure  15:  Variation  of  skewness  with  vertical  distance  from  shear  layer. 


c 


2 


7 P 
P 


(19) 


Using  conservation  of  mass  to  replace  the  material  derivative  of  density,  the  equation  for  pressure  is 


Dp 

~Dt 


dui 


p  DS 
cv  Dt 


acoustic  entropic 


(20) 


For  now  we  will  focus  on  the  acoustic  contributions  to  the  pressure  changes.  To  assess  the  magnitude 
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Figure  16:  Fluctuation  pressure  trace  at  =  5.5  for  (a)  M  =  1.5,  (b)  2.5,  and  (c)  3.5. 


x  x 

Figure  17:  (a)  Schematic  of  the  very  near  field  acoustic  field  with  V  •  u  <  —0.0125  below  the  vorticity 
magnitude  colored  between  |V  x  u|  =  0.1  and  5A U/5m(t).  (b)  Detailed  view  of  wavefront  with  normal 
vectors  shown  at  each  computational  grid  point. 


of  the  entropic  contributions,  the  left-hand-side  of  equation  20  will  be  directly  computed  from  our  DNS 
database.  To  find  an  evolution  equation  for  pressure  variance  ((p')2)  and  skewness  (( p ')3)>  the  variables 
must  be  Reynolds  decomposed  by 


p  =  p'+p  (21) 

Ui  =  u[  +  Wi , 

where  ()'  and  ()  denote  a  fluctuation  and  an  average.  The  method  of  averaging  depends  on  the  problem. 
For  a  parallel  shear  flow,  the  average  is  taken  to  be  an  ensemble  average  over  the  x-z  plane  at  a  given  y. 
Substituting  the  Reynolds  decomposition,  the  pressure  equation  is 


Dp +  p' 
Dt 


-7  (p  +  P') 


dui  +  v!i 
dxi 


+  . . . , 


(22) 
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Figure  18:  M  =  1.5:  (a)  Dilatation,  (V •  u)5m(t) / AU <—0.0125,  when  <5m(t)/<5^  =  32.8  at  z  =  Lz/2.  (b) 
Average  wave  angles  ((d))  conditioned  on  (V  •  u)Sm(t) / AU <—0.0125  shaded  by  one  standard  deviation 
(a)  about  the  mean. 
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Figure  19:  M  =  2.5:  (a)  Dilatation,  (V  ■  u)<Sm(i)/A!7 <—0.0125,  when  (5m(t)/<5^  =  32.8  at  2  =  Lz/ 2.  (b) 
Average  wave  angles  ((d))  conditioned  on  (V  •  u)Sm(t)  /  AU <—0.0125  shaded  by  one  standard  deviation 
(cr)  about  the  mean. 


where  (. . .)  represent  the  entropic  contribution  which  will  not  be  explicitly  written  out  here  but  will  com¬ 
puted  indirectly  by  evaluating  the  left-hand-side  from  the  DNS  database.  For  the  rest  of  the  exposition, 
the  entropic  components  will  be  left  out  of  the  equations.  Taking  the  average  equation  22,  the  average 
pressure  evolution  is 


dp  _ dp 

dt  1  dxi 


,  dp' 

“•as;"-7 


dui 


P-. 


-P 


,duj 


dxi  dxi 

Again  using  Reynolds  decomposition,  the  evolution  of  fluctuating  pressure  is  given  as 

dp'  dp  dp 
~dt  ~  ~dt~  ~dt' 

and  by  substituting  the  full  and  mean  pressure  equation,  the  fluctuating  pressure  is 


dp'  _ dp'  ,  dp  ,  dp'  dp' 

Ik  =  ~Uldx)  ~Uldx,~Uldxj.+Uldx(i. 


■  7 


_du(  ,dui  ,  du't 


dxi 


-P 


dxi 


-P 


dxi 


+  7 p' 


duj 
dx  1 


(23) 


(24) 


(25) 


With  the  evolution  equation  for  fluctuation  pressure,  the  high-order  moment  equations  for  the  variance 
skewness  can  be  written  as 
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Figure  20:  M  =  3.5:  (a)  Dilatation,  (V •  u)5m(t) / AU <—0.0125,  when  =  32.8  at  z  =  Lz/2.  (b) 

Average  wave  angles  ((d))  conditioned  on  (V  •  u)Sm(t) / AU <—0.0125  shaded  by  one  standard  deviation 
(a)  about  the  mean. 


Figure  21:  Variation  of  average  wave  angle  with  dilatation  threshold  for  (a)  M  =  1.5,  (b)  2.5,  and  (c) 
3.5.  The  various  values  of  dc  shown  range  between  -0.0125  and  -0.003125. 
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Figure  22:  Scaled  average  number  of  waves  conditioned  on  (V  •  \i)8m(t)  /  AU  <  —0.0125. 


d(P')2  _  ,dv' 

-dT  =  v^i 


d(p')3 

dt 


(26) 
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Figure  23:  Variation  of  scaled  average  number  of  waves  with  dilatation  threshold  for  (a)  M  =  1.5,  (b) 
2.5,  and  (c)  3.5.  The  various  values  of  dc  shown  range  between  -0.0125  and  -0.003125. 


respectively.  Multiplying  the  equation  25  by  2 p'  and  taking  the  average,  the  pressure  variance  equation 
is 


d(p')2  0_  ,  dp'  0  ,  ,  dp' 

wr  =  -2 Uipf- - 2  u'j/  -  2w'y — 

at  oxj.  oxi  oxi 


-27 


_dm  dm —-2  ,  A29u- 

+  &T(p  )  +(P) 

L/  J/  7  L/  J/  o  7 


Moving  the  first  term  on  the  right  hand  side  to  the  left  we  can  define 


d{pl)2  I  -d{v'f 
Dt  dt  1  dxi 


leaving 


D -  dp -  dp' 

-(pT  =  -2g£«ip'-2«yf|--27 


_dm  duiT—^  du) 

pd7p +  d^{p) +  {p) d^ 

L'il'7  1/ 7  L/1I/7 


(27) 


(28) 


The  equation  for  pressure  variance  differs  from  Sarkar  (1992)  (see  equation  7)  and  could  be  made  identical 
by  invoking  the  following  relation 


^  =  2 M?f  +  wdA. 

dxi  dxi  dxi 

Performing  the  same  manipulations  for  pressure  skewness  we  have 


(29) 


1  dip')3  1  d(p')3 

3Adr  +  r^  =  -^u: 


Using  a  simplifying  relation 


the  equation  for  skewness  is 


,du' 


.du'r 


-^y)2  -  £(p,)3  -  (p,)ag +p,^y)2 


dp'p'p'K  _  2  ,dp^ 

Ox,  ~6{P)  U'dXi  +[P)  dxi’ 


(30) 


(31) 


D  3  q  dP  !(  A2  dip'fu) 

DtV)  =-^uW) 


,  du'i 


,  dp' 


dxi 


+  (J»03g7+3(!>')%' 


+37 


dui 


dui 


,du' 


.du'i 


-p-x—{p')2  -  ^—{p')3  -  ( P 03^  +P'xxA{p')2 
dxi  dxi  dxi  dxi 


(32) 
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Figure  24:  Visualizations  of  merging  wave  fronts  for  M  =  2.5  when  Srn (t) / &°n  =  15  to  25  at  z  =  Lz/ 2: 
color  contours  show  |V  x  u|  between  0.1  and  5A U/5m(t)  and  gray-scale  contours  show  V  •  u  between 
±0.5AU/5m{t). 


- 3/2 

The  entire  equation  can  be  normalized  by  ( p ')2  to  give  the  evolution  of  normalized  pressure  skewness, 
Skip1)-  The  skewness  and  variance  evolution  equations  are  general  for  any  flow  configuration.  We  will 
specialize  them  to  a  parallel  shear  flow  where  only  derivatives  of  mean  quantities  in  the  y-direction  are 
nonzero. 


DT7vi  odP  77  02  dip'YV 

mw)  =~Ww)  ST 


+  ( p’)3d '  +  3  {p'Yui 


,  dp' 


dxj. 


+3q 


dv- 


—pd'ip')2  —  t^-(p')3  —  ( p')3d '  +p'd'(p')2 
dy 


where  d'  =  dv!ijdxi  represents  the  fluctuating  dilatation. 


(33) 
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Figure  25:  (a)  Schematic  of  space-space  correlation,  (b)  Integral  length  when  =  35. 
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Figure  26:  Space-space  correlation  contours  of  pressure  at  y/5m{t)  =  15  when  =  35  for  (a) 

M  =  1.5,  (b)  2.5,  and  (c)  3.5.  The  dashed  line  corresponds  to  the  orientation  with  respect  the  horizontal 
axis  of  the  maximum  integral  length  scale  of  the  two-dimensional  correlation  function. 


5.2  Results 

For  the  DNS  of  M  =  0.9  and  M  =  2.5,  the  skewness  and  variance  budgets  of  pressure  are  computed. 
First,  it  is  important  to  assess  our  approximation  of  the  time  rate-of-change  of  variance  and  skewness. 
A  first-order  approximation  of  the  derivative  of  variance  is  given  by 

dip')2  _  (p')2(t  +  At)-  jp')2it) 

dt  ~  At  1  ’ 

where  the  pressure  field  p'{t  +  At)  is  evaluated  using  the  standard  fourth-order  Runge-Kutta  by  solving 
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the  entire  compressible  Navier-Stokes  equations  at  a  future  time  of  t  +  At.  The  same  approach  is  used  for 
the  skewness.  Figure  27  shows  the  first-order  error  of  approximating  the  variance  and  skewness  temporal 
derivatives  at  the  centerline  (y  =  0) .  As  the  time  step  is  reduced,  the  approximation  follows  the  expected 
—  1  slope  but  then  diverges  due  to  round-off  errors  associated  with  evaluating  the  full  equations  at  time 
steps  that  are  numerically  very  close.  A  At  =  10~5  is  used  to  approximate  the  rate  of  change  of  variance 
and  skewness. 


Figure  27:  Error  in  first-order  approximation  of  time  rate-of-cliange  of  variance  and  skewness. 

Figures  28  and  29  show  the  pressure  variance  budget  for  different  y  locations  in  the  near  acoustic 
field  as  a  function  of  time  for  M  =  2.5  and  0.9,  respectively.  The  statistics  were  averaged  at  the  ±y. 
The  time  traces  of  the  components  of  the  budget  show  how  the  initial  transitioning  shear  layer  radiates 
sound  for  t  <  100.  For  both  cases,  the  pressure-diliation  term  makes  up  most  of  the  acoustic  distribution. 
For  the  M  =  2.5  case,  the  entropic  contribution  is  slightly  larger  than  the  M  =  0.9  which  is  expected 
since  at  M  =  2.5  contain  Mach-like  waves  which  are  weak  shock-like  features  with  steep  gradients  that 
can  dissipate  energy.  For  each  y  location  the  time  rate-of-change  of  variance  in  the  time  average  sense 
is  larger  than  zero.  This  suggests  that  the  intensity  grows  with  time  which  is  expected  to  occur  for  a 
growing  shear  layer  whose  volume  increases  with  time  meaning  it  will  radiate  more  sound. 

Figures  30  and  31  show  the  pressure  skewness  budget  for  different  y  locations  in  the  near  acoustic 
field  as  a  function  of  time  for  M  =  2.5  and  0.9,  respectively.  For  each  y  location  the  time  rate-of-change 
of  skewness  in  the  time  average  sense  is  nearly  zero.  This  suggests  that  as  the  shear  layer  grows  the 
normalized  skewness  will  remain  constant. 


6  Conclusions 

DNS  of  high-speed  turbulent  shear  layers  for  M  =  0.9,  1.5,  2.5,  and  3.5  have  been  completed.  Analysis 
of  acoustics  show  near-field  nonlinearity  of  wave  merging  and  flattening  in  the  near  acoustic  field  are 
important  for  skewness.  Analysis  of  the  skewness  budget  suggests  that  acoustic  interaction  effects  are 
balanced  by  entropic  contributions.  This  means  the  initial  footprint  of  skewness  must  arise  near  the 
turbulence  source.  Skewness  variation  with  distance  from  the  shear  layer  supports  this  claim. 

7  Personnel  supported  during  duration  of  grant 

The  work  is  being  supervised  by  Prof.  Freund  at  UIUC.  Funding  supported  the  PhD  studies  Mr.  Aaron 
Anderson  and  Mr.  David  Buchta. 


26 


Figure  28:  Budget  of  pressure  variance  for  M  =  2.5  at  y/5m(t)  =  7.5  (a),  10  (b),  and  15  (c).  All 
components  have  been  normalized  by  {p')2  ■ 

8  AFOSR  point  of  contact 
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A  Flow  equations 

The  nondimensional  indpendent  variables  of  space  and  time  are 


x  y 

x  =  — ,  y  = 


l* 


l*'  Z  l* 


t  = 


t*c* 


and  the  flow  variables  have  been  nondimensionalized  by 


(35) 
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Figure  31:  Budget  of  pressure  skewness  for  M  =  0.9  at  y/Sm(t )  =  7.5  (a),  10  (b),  and  15  (c).All 
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components  have  been  normalized  by  ( p ')2 


n*  ii*  T  c 

P  Ui  rri  ±  cp 

P  =  - ,  Ui  =  — ,  1  =  t - 5r,  p  =  — 

o*  c*  ( c*  o* 

roo  °oo  v'-'oo  /  ro 


p' 

fc*  l2' 

•  \L,ooJ 


(36) 


The  ()*  quantities  represent  dimensional  ones  and  the  freestream  values  of  density,  pressure,  and  speed 
of  sound  are  denoted  by  p^ ,  Poo ,  and  Coo,  respectively.  The  specific  heat  at  constant  pressure  is  cp. 
With  the  nondimensionalization  given  above,  the  compressible  flow  equations  of  mass,  momentum,  and 
energy  are 


dp  d(puj) 
dt  dxi 


(37) 
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(38) 


de 

dt 


d(pui) 

dt 


_d_ 

dxi 


( PUiUj ) 


dp 

dxi 
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dxi 


[ui(e  +  p)] 


7  d2 
RePr( 7  —  1)  dxf 


drij_ 
dxi  ' 


dxi 


{pj  Uj  ) , 


(39) 


respectively.  The  viscous  stress  tensor  assuming  a  Newtonian  fluid  and  neglecting  the  effect  of  bulk 
viscosity  is  expressed 


_  1  /  dm 

Re  \  dxj 

The  total  energy  assuming  an  ideal  gas  law  is 


dui 

dxi 


2  duk  \ 

3  dxk)  ’ 


(40) 


P  1 

e  = - -  +  -pu.iUi. 

7  —  1  2 


(41) 


The  flow  equations  are  parameterized  by  the  ratio  of  specific  heats,  Reynolds  number,  Prandtl  number 
given  by 


7  =  %  =  1-4  (42) 

^ v 

n*  r*  7* 

Re  =  Po°  °°  (43) 

p* 

ctu* 

Pr  =  —  =  0.7.  (44) 

K* 

The  dynamic  viscosity  (p*)  and  thermal  conductivity  (k*)  are  assumed  constant.  The  Reynolds  number 
can  be  rescaled  to  reflect  a  more  useful  comparison  between  flows  of  different  speed.  The  Reynolds 
number  based  on  the  momentum  thickness  and  velocity  difference  is  written  as 


R-esm 


PooAU*S*m(t) 

p* 


(45) 


The  relationship  between  the  simulation  Reynolds  number  and  momentum  thickness  Reynolds  number 
is 


where  the  Mach  number  is 


ReSm  =  ReM6*m{t)/l*, 


M  = 


(46) 


(47) 


B  Numerics 


B.l  Discrete  Fourier  method 

The  discrete  Fourier  transform  of  a  periodic  function  /  represented  by  N  uniformly  spaced  grid  points  is 


f(kj)  =  —  f(xm)  exp(— v7— Ifcja;TO), 

m—0 

where  xm  is  the  spatial  location  at  the  mth-grid  point  and  the  wavenumber  is  defined  as 

2irl 

1  ~  NAx' 

The  discrete  inverse  Fourier  transform  is  given  as 


(48) 


(49) 
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(50) 


N/  2  —  1 

f(xj)  =  f(krn)exp{V^ikmXj), 

m—N/2 

and  the  first  and  second  order  derivatives  can  be  written  as 

N/2-1 

fix,)  =  ^2  V^ikmf(km)exp(V^ikmXj)  and  (51) 

m=N/2 

N/  2-1 

f"{xj)  =  X]  _fcm/(fcm)exp(v/zIfcTOXj),  (52) 

m—N/2 

respectively. 


B.2  Finite  difference  first  derivative 

The  approximation  of  the  first  derivative  at  spatial  location  xQ  can  be  written  as 
d  M  i  N 

-q^J(xo)~  X!  bif'(xo+iAx)  =  —  ^2  ajfixo+j Ax),  (53) 

i=—M  j——N 

where  the  coefficients  are  assumed  to  have  symmetry  6,;  =  b—i  and  a*  =  —  a_j.  Traditionally  (Lele, 
1992),  the  coefficients  are  determined  by  matching  terms  in  the  Taylor  series  substituted  in  equation  53. 
Others  (Tam  &  Webb,  1993;  Bogey  &  Bailly,  2004)  have  used  a  combined  approach  to  determine  the 
coefficients.  The  truncated  Taylor  series  approximation  is  used  to  enforce  a  desired  order  of  accuracy,  and 
the  remaining  coefficients  are  determined  by  minimizing  the  error  between  the  exact  and  approximate 
dispersion  relation.  The  properties  of  the  scheme  then  preserve  aspects  of  the  dispersion  properties  of 
the  derivative.  For  an  explicit  (M  =  0)  7-point  stencil  ( N  =  3)  each  of  the  terms  in  finite  difference 
expression  is  written  as  a  Taylor  series  expansion  up  to  0(Ax6)  by 
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Subsitituting  these  relations  into  the  differencing  equation  and  matching  like  terms,  the  constraint  equa¬ 
tions  that  enforce  formal  order  of  accuracy  up  to  sixth  order  are 


a0 

=  0 

(54) 

2a\  +  4a2  +  6a3 

=  1  (second  order) 

(55) 

2 

2  •  23  2  •  33 

3!ai  + 

3!  “2  +  3!  “3 

=  0  (fourth  order) 

(56) 

2 

2  •  25  2  ■  35 

5!ai  + 

5!  02  +  5!  “3 

=  0  (sixth  order). 

(57) 
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Directly  solving  of  equations  55-57  results  the  standard  sixth-order  approximation  with  coefficients 


a0  =  0 

at  =  -a- 1  =  0.750000000000000 

a2  =  -a- 2  =  -.150000000000000 

a3  =  —  a_3  =  0.0166666666666667. 


Up  to  this  point,  there  has  been  no  guarantee  that  the  finite-difference  approximation  mimics  the  prop¬ 
erties  of  the  first  derivative  in  wavenumber  space.  To  obtain  a  resolution-optimized  fourth-order  method 
for  example,  relations  55  and  56  are  enforced  leaving  one  free-parameter  to  be  tuned.  The  tuning  scheme 
involves  an  optimization  in  Fourier  space.  The  scheme  can  be  written  in  Fourier  space  by  substitut¬ 
ing  the  inverse  discrete  Fourier  transform  (equation  50).  Ignoring  the  summation  over  the  modes  the 
differencing  approximation  is 

M  N 

bi\/—lkf(k)  exp  [V— 1  k(x  +  *Ax)]  =  ajf(k)  exp  [y/—lk(x  +  jAx)]  ,  (58) 

i——M  X  j=—N 

and  cancelling  f(k)  exp  [v/— 1 k(x)]  from  both  sides  we  have 


M 

biV—  lfcexp  \V—  lfc(*Aa;)] 

i—  —  M 

Isolating  k Ax  on  the  left-hand-side 


1  J< ' 

—  ^  a,j  exp  [V^lk(jAx)]  . 

j——N 


kAx 


and  recalling  Euler’s  formula 


N 

-v73!  E  ajexp[y/^lk(jAx)\ 

j=-N 

M 

E  h  exp  [v/— Tfc(iAa:)] 

i——M 


(59) 


(60) 


cos(x) 
sin  (a; ) 

the  effective  wavenumber  equation  is 


exp(\/— la;)  +  exp(—  \/—  lx 


2 

exp(v/— la:)  —  exp(—  y/—  lx 

2v/=I 


N 

2  V  a,-  sinOfcAa;) 

(kAx)'  -  - .  (61) 

b0  +  2  E  bi  cos(ikAx ) 

i= 1 

The  explicit  finite-difference  (M  =  0)  eciuation  approximates  the  first  derivative  as  a  truncated  sine 
series  in  Fourier  space.  A  real-valued  expression  for  equation  61  is  a  consequence  that  the  coefficients 
were  symmetric;  otherwise  the  dispersion  relation  would  be  complex.  The  (kAx)'  notation  in  equation 
61  denotes  an  approximation  to  the  exact  wavenumber  relation  (kAx). 

To  preserve  aspects  of  the  actual  dispersion  relation,  the  integral  error, 


-ln(6) 


E  = 


[(kAx)  —  (kAx)']  dln(kAx), 


’  ln(a) 


(62) 


between  the  exact  and  effective  wavenumber  will  be  minimized  to  obtain  a  final  relation  to  solve  the  final 
tuning  parameter.  The  limits  of  integration  for  a  seven-point  M  =  3  scheme  is  a  =  7t/16  and  b  =  7r/2. 
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The  thirteen-point  schemes  uses  a  =  7r / 16  and  b  =  3tt/5.  For  the  fourth-order  M  =  3  scheme,  the  final 
constraint  that  optimizes  the  integral  error  is 


dE 

da3 


=  0. 


(63) 


Figure  33  shows  the  effective  wavenumber  and  associated  error  with  respect  to  the  exact  relation.  By 
sacrificing  order  of  accuracy  for  resolution,  the  optimized  formulas  (abbreviated  opt.)  have  more  error 
in  the  low-wavenumber  range.  The  standard  schemes  (abbreviated  std.)  have  coefficients  that  maximize 
formal  order  of  accuracy.  For  the  seven-point  M  =  3,  fourth-order  resolution  optimized  scheme  the 
coefficients  are 


a0  =  0 

ai  =  — a_i  =  0.824639848382100 
a2  =  — a_2  =  —0.209711878705680 
a3  =  — a_3  =  0.0315946363430865 

The  coefficients  for  the  M  =  6,  fourth-order  scheme  are 


a0  =  0 

CL\  =  — CL— i 

Cl2  =  — a_2 

a3  =  -a_  3 

ai  =  — a_4 

a5  =  — a_5 

a6  =  -a- 6 


0.9108405208695360 

-0.3419538377082619 

0.1380399894807369 

-0.04827039810294327 

0.008624243759248696 

-0.0006681390707308727. 


Figure  32:  Effective  wavenumber  of  the  interior  first  derivative  (a)  and  associated  difference  from  the 
exact  (b). 

Near-boundary  points  use  reduced-order  central  differences  as  described  in  Freund  (1997a).  On  the 
boundaries  a  third-order  one-sided  approximation  is  used. 


B.3  Finite  difference  second  derivative 

The  explicit  finite  difference  approximation  of  the  second  derivative  is  given  by 

(92  l  N 

g^f(xo)  ~  f"(x0)  =  Y,  cjfixo  +  jAx),  (64) 
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with  constraint  equations  for  N  =  3  on  coefficients  enforcing  formal  order  of  accuracy  up  to  sixth-order 

by 


Directly  solving  equations 


Co  H-  2ci  H-  2c2  2c3 

=  0 

(65) 

2 

2  •  22  2  •  32 

=  1  (second  order) 

(66) 

2!Cl 

+  2!  c2+  2,  C3 

2 

2  •  24  2  •  34 

=  0  (fourth  order) 

(67) 

if' 

+  4!  C2+  4!  C3 

2 

2  •  26  2  •  36 

=  0  (sixth  order). 

(68) 

6!Cl 

+  6!  C2+  6!  C3 

(69) 

65-68 

results  in  the  standard  sixth-order  scheme  with  coefficients 

c0  =  —2.72222222222222 

ci  =  c_i  =  1.50000000000000 
c2  =  c_2  =  -.150000000000000 
C3  =  C_3  =0.0111111111111111. 

Substituting  the  discrete  Fourier  representation  of  the  function  in  equation  72,  a  modified  wavenumber 
expression  for  the  explicit  finite-difference  approximation  of  the  second  derivative  is 


N 

( kAx )"  =  ca  +  2  Cj  cos(jAx).  (70) 

i=i 

The  modified  wave  number  approximates  the  exact  dispersion  relation  by  a  truncated  cosine  series.  The 
square  error  between  the  exact  and  effective  wavenumber  on  a  scaled  periodic  domain  [0,  n)  is  given  by 

pln(b) 

E=  I  [(kAx)2  —  (kAx)"] 2  dln(kAx),  (71) 

J  ln(a ) 

where  the  bounds  of  integration  for  a  seven-point  scheme  are  a  =  7r / 16  and  b  =  7t/2.  For  a  thirteen-point 
scheme  a  =  ir/16  and  b  =  3tt/5.  Figure  33  gives  the  modified  wave  number  and  error  associated  with 
seven-point  and  thirteen-point  schemes  of  the  second  derivative.  Like  the  first  derivative,  the  optimized 
methods,  by  sacrificing  maximum  accuracy  for  resolution,  lower  wavenumbers  are  increased  error. 


(b) 


Figure  33:  Effective  wavenumber  of  the  interior  second  derivative  (a)  and  associated  difference  from  the 
exact  (b). 


The  coefficients  for  the  resolution  optimized  M  =  3  and  6  schemes  are 
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Co  =  -2.81321312844389 
d  =  c_  i  =  1.56824317966625 
c2  =  c2  =  -.177297271866500 
c3  =  c3  —  0.0156606564221945. 


c0  =  -3.093268941798102 
ci  =  c_i  =  1.811759211319885 
c2  =  c2  —  -.3343827021872651 
c3  =  c3  =  0.08734940331848545 
c4  =  c4  =  -0.02186933368500581 
c5  =  c5  =  0.004224346387577600 
c6  =  c6  =  -0.0004464542546261400. 


respectively.  Near-boundary  points  use  reduced-order  central  differences.  For  i 
standard  fourth-order  method  is  used  with  coefficients 


2,  *  =  N  —  2,  a 


Co  —  —2 

ci  =  a_i  =  1 


For  i  =  1,  i  =  N  —  1,  second-order  scheme  is  used  with  coefficients 


co  = 

4 

Cl  =  C_  1  =  3 


On  the  boundaries,  a  one-sided,  third-order  method  is  used  where  the  coefficents  at  i  =  1  are  given  by 


Co 

=  2 

Cl 

=  -5 

C2 

=  4 

C2 

=  -1 

B.4  Selective  high-wave  number  filtering 

It  may  be  necessary  to  apply  a  selective  higli-wave  number  filtering  to  the  solution  so  that  parasitic, 
short  wavelength  waves  do  not  deteroriate  the  solution.  The  differencing  equation  is  similar  to  Bogey  & 
Bailly  (2004)  given  here  as 

N 

f{x0)  =  f(x0)-  X!  hjf(xo+jAx).  (72) 

j=-N 

For  a  thirteen-point  (N  =  6)  scheme  and  assuming  hi  =  h- 1,  the  constraint  equations  of  formal  order 
of  accuracy  are  written 
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A  transfer  function  in  wavenumber  space  can  be  found,  much  like  the  modified  wavenumbers  of  the 
derivative  approximations.  For  the  explicit  formuals  used  here,  the  transfer  function  is 


T(kA)  =  1  —  +  2  ^  hj  cos(jAa;)  J  , 

and  the  square  integral  error  between  the  filtered  and  un-filtered  solutions  is 

rln(h) 


E  = 


[1  —  T(kAx)\  dln(kAx), 


(80) 


(81) 


'  ln(a 


where  the  limits  of  integration  for  the  seven-point  are  a  =  7r/16  and  b  =  37r/32.  For  the  thirteen-point 
formulas  the  limits  are  a  =  7t/16  and  b  =  7r/2.  The  coefficients  for  the  seven-  and  thirteen-point  formulas 
by  minimizing  equation  81  are  given  as 


h0  =  0.288090709024785 
hi  =  h- 1  =  -.228272677256196 
h2  =  h_ 2  =  0.105954645487608 
h3  =  h- 3  =  -0.02172732274380381 


h0  =  0.1913963126291382 
hi  =  h_x  =  -0.1719104977994496 
h2  =  h—2  =  0.1237396648728178 
h3  =  h- 3  =  -0.06976715092480617 
hi  =  h—i  =  0.02936617293645158 
h5  =  h-  5  =  -0.008322351275744250 
h6  =  h—o  =  0.001196005876161525 


respectively. 

When  an  implicit  filtering  scheme  is  desired,  the  finite  difference  equation  can  be  written  as 

M  N 

^  gif(x0  +  iAx)=  y  hjf(x0+  jAx),  (82) 

i——M  j——N 
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E-l 


(b) 


Figure  34:  Transfer  function  of  the  interior  high-wavenumber  filter  (a)  and  associated  difference  from 
the  exact  (b). 


and  the  same  process  as  described  above  for  selecting  the  stencil  width,  formal  order  of  accuracy,  and 
enforcing  an  effective  wavenumber  behavior  is  carried  out.  Following  Lele  (1992),  the  resolution  behavior 
in  wavenumber  space  is  set  differently  than  Tam  &  Webb  (1993);  Bogey  &  Bailly  (2004).  Values  of 
T(fcA)  at  a  particular  fcA  to  find  the  remaining  degrees  of  freedom  of  the  stencil  formula.  The  behavior 
at  fcA  =  7r  is  also  set  to 


G?2T(fcA) 

d2kA 


(tt)  =  0. 


(83) 


Though  this  approach  does  not  minimize  the  integral  square  difference  between  the  actual  and  effective 
wave  number  resolution,  the  resolution  properties  are  well-known  and  have  not  had  any  negative  effects 
on  the  solution.  The  trends  for  the  seven-point  implicit  scheme  (C.2.10.b)  designed  by  Lele  (1992)  are 
given  in  figure  34. 


C  Reynolds  number  effect 

A  higher  Reynolds  number  simulation  at  Re-s^  =  120  was  run  for  M  =  2.5.  The  grid  in  the  x-,  y-  and 
z-directions  was  uniformly  spaced  with  Ax  =  Ay  =  A z  =  0.5(5)^  which  is  nearly  two  times  the  resolution 
of  the  Res^  =  60  simulations.  The  numerics  of  the  high-Reynolds  number  simulation  were  modified.  The 
first-order  derivatives  in  all  directions  were  computed  with  a  13-point,  fourth-order  resolution  optimized 
scheme  designed  by  Bogey  &  Bailly  (2004).  A  13-point  fourth-order  filter  was  applied  every  5  time 
steps  in  all  three  coordinate  directions.  The  second  derivative  was  evaluated  by  repeated  first-order 
derivatives.  High-order,  non-centered  schemes  crafted  by  Berland  et  al.  (2007)  were  used  for  derivative 
approximations  near  boundaries. 

An  identical  numerical  scheme  in  all  coordinate  directions  of  equal  spatial  resolution  tests  the  robust¬ 
ness  of  the  initial  numerical  scheme.  This  verifies  the  sufficient  resolution  of  the  initial  simulations  and 
that  any  potentially  negative  anisotropic  effects  of  the  numerical  method  did  not  affect  the  statistics  of 
interest. 


D  Instability  modes  of  inviscid  parallel  temporal  shear  layers 

The  primative  flow  variables  are  expressed  by  a  Reynolds  decomposition, 

[u,  v,  w,  p,  T }  =  [u,  v,  w,  p,  T]  +  [u ',  v',w',  p\  T'] ,  (84) 

where  ()  is  an  average  in  the  statistically  homogeneous  directions  and  (')  represent  a  perturbation  quan¬ 
tity.  The  mean  flow  quantities  are  prescribed  by 
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[u,v,w,p,T](y)  =  [u(y),  0, 0, 1  ,T(u)] , 
with  small-amplitude  (e)  perturbations  in  the  form  of  travelling  waves 


(85) 


[u',v',w',p',r] 


eRe  I 


u,  v,  w,  p 


T 


(■ y )  exp{i(aa;  +  ftz  —  . 


(86) 


This  same  formulation  has  been  described  elsewhere  (Sandham,  1989)  and  a  brief  outline  will  be  explained 
here.  Substituting  the  mean  and  perturbation  quantities  into  the  invsicid,  compressible  flow  equations 
and  neglecting  quadrating  interactions  of  perturbation  quantities,  a  set  of  ordinary  differential  equations 
is  given  as 


,  „ dp 

piyau  —  oj)  +  v— — b  p 

ay 


. ,  „  n  -  \  dv 

i  (au  +  pw)  +  — 

ay 

. ,  _  ,  .  ,  du 

i  (au  —  oj)u  +  v  — 

ay 

p  [i(au  —  u>)v] 
p  [i(aw  —  ui)w\ 


i  (au  —  uj)T  + 


(IT 
dy 

-(7-1) 


P 


=  0 

— lap 
=  ^M? 
dp  1 

"  dy-U't 
=  ~i  ftp 
7  Ml 


i  (au  +  0w) 
=  pT  +  pT 


dv 

dy 
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(88) 

(89) 

(90) 

(91) 

(92) 


Using  the  method  of  Gropengiesser  (1970),  the  equations  can  be  combined  to  form  two  equations  for 
v  and  p  and  then  combined  into  a  single  eigenfunction  in  the  from 


X  = 


lap 


with  evolution  equation 


where 


pMlv  ’ 

dx  a2(u  —  u;/a)  x(X9  +  ^) 


dy 


T 


(u  —  ui/a) 


a2  +  /32  M2(au  —  ui )2 


9=  -  2  2 

paz  az 

By  inspecting  the  behavior  of  the  v  and  p  eigenfunctions  at  y  — >  ±oo,  the  behavior  of  \  is 

a(u  —  ui/a) 


x(y  =  ±oo)  = 


VgT 


(93) 


(94) 


(95) 


(96) 


The  mean  velocity  profile  as  set  as 


AU  , 
u  =  — —  tanh 


(97) 


and  the  temperature  relation  is  specified  from  the  Crocco-Busemann  relation  assuming  unity  Prandtl 
number.  For  temporally  developing  shear  layers,  a  real-valued  mode  (a,  (3)  is  specified  with  an  initial 
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Figure  35:  Variation  of  two-dimensional  unstable  modes  with  streamwise  wavenumber,  a:  (a)  temporal 
amplification  rate  and  (b)  rate  of  phase  change. 


guess  of  complex- valued  ui.  The  evolution  equation  for  \  in  equation  (94)  is  integrated  using  a  fifth-order 
adaptive-step-refinement  Runge-Kutta  scheme  (Press  et  al.,  1986)  from  free-stream  boundary  conditions 
in  equation  (96)  to  the  centerline  at  y  =  0.  The  solutions  at  y  =  0  from  each  of  the  free-stream  conditions 
are  compared  and  an  update  to  the  initial  ui  guess  is  made  using  a  complex-valued  Newton-Raphson 
method.  Once  the  two  values  of  \  match  at  y  =  0  to  within  specified  tolerance  of  10”14,  the  procedure 
is  stopped. 

Figure  35  shows  the  solution  of  ui  for  two-dimensional,  spanwise  oriented  modes  (a,  0).  As  the  Mach 
number  increases,  the  amplification  rate  decreases.  For  Mi  >  2,  the  phase  speed  of  the  mode  is  finite. 
These  two-dimensional  modes  are  known  as  the  supersonic  modes  (Sandham,  1989),  moving  supersonic 
with  respect  to  the  free-stream  velocity. 

The  behavior  of  three-dimensional  amplification  rate  are  shown  in  figure  36.  As  the  Mach  number 
increases,  the  amplification  rate  decreases.  For  M  >  1.5,  the  most  amplified  mode  is  three-dimensional, 
whereas  for  M  <  1,  the  two-dimensional  modes  are  the  most  amplified. 

For  M  =  0.9, 1.5,  2.5,  3.5,  the  angle 

0  =  tan-1  ) 

corresponding  to  the  most  unstable  mode  is  given  in  figure 
suggestively  follows  the  relation 

Mc  cos(0)  =  0.6,  (99) 

an  observation  from  the  stability  calculations  of  Sandham  &  Reynolds  (1990).  The  collapse  of  this  data 
with  Mc  is  rather  robust  to  effects  of  varying  velocity  ratio,  temperature  ratio,  and  free-stream  stagnation 
entalpies  (Sandham  &  Reynolds,  1990).  For  Mc  >  0.6,  the  Mach  number  normal  to  the  unstable  wave 


(98) 

37.  For  Mc  >  0.6,  the  trend  of  the  angle 
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0  =  tan  1(/3/a) 


Figure  36:  Variation  of  three-dimensional  growth  rates  with  angle  between  the  streamwise  and  spanwise 
wavenumbers:  (a)  a  =  0.25  and  (b)  0.5. 


Figure  37:  Angle  of  the  most  unstable  three-dimensional  mode  with  Mc.  The  line  is  an  empirical  relation 
given  by  Sandham  &  Reynolds  (1990)  for  Mc  >  0.6. 
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crests  is  M  =  0.6.  From  the  perspective  of  the  sonic  eddy  model  of  Breidenthal  (1990)  to  participate  in 
energy  transfer  processes,  the  subsonic  speed  normal  to  these  waves  that  is  predicted  by  linear  theory 
suggests  this  is  an  important  mechanism  for  instability  of  these  shear  layers  at  high  speeds. 

A  connection  can  be  made  between  the  orientation  of  the  wavefronts  and  the  theoretical  linear  phase 
speed.  The  wavelength  between  wave  crests  in  the  y-direction  (ly)  and  the  streamwize  wavelength  of  the 
unstable  mode  ( lx  =  2ir/kx)  (shown  in  figure  D)  provides  a  measure  of  the  orientation  of  the  waves  given 

by 


0  =  tan  1 


=  tan  1 


/  17.87  \ 
V2tt/0.289  ) 


39.42°. 


(100) 


Figure  38:  Schematic  of  radiating  pressure  eigenfunction  for  M  =  2.5  with  kx  =  0.289  and  ly  =  17.87. 


This  direct  measurement  is  compared  to  an  inferred  orientation  based  on  the  phase  speed  of  the 
mode  relative  to  the  free-stream  velocity.  The  phase  speed  is  expressed  as  cm  =  Re(uj)/a  =  0.259. 
Given  the  non-dimensionalization  of  the  disturbance  equations  this  must  be  re-scaled  by  the  free-stream 
Mach  number  M  =  1.25  giving  a  velocity  of  0.323.  The  wave  angle  is  then 

em  =  sin  1  (^  — )  =  sin  1  (ym  +  v-J  =  39.46  .  (101) 

The  agreement  between  the  two  approaches  is  promising;  though,  linear  analysis  overpredicts  the 
wavefront  orientation.  For  M  =  2.5,  DNS  predicted  an  average  angle  of  6  =  34.2°  at  5TO(t)/^  =  15, 
suggesting  structures  are  convecting  faster  than  predicted  by  linear  analysis.  Empirical  relations  from 
Oertel  (1979)  suggests  dominance  of  waves  at  an  orientation  of 

6  =  sin-1  (t^T2j)  =  34-85°>  (102) 

assuming  the  equal  free-stream  sound  speeds.  This  relation  seems  to  correspond  more  closely  to  the 
observed  orientation  in  the  DNS  data  which  may  not  be  surprising  since  the  Oertel  (1979)  relations  are 
correlated  to  turbulent  supersonic  shear  layers.  Simulation  by  Nichols  et  al.  (2013)  have  also  shown 
qualitative  agreement  of  wavefront  orientation  with  using  this  relationship. 
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Abstract 

Fighter  jets  and  other  aircraft  with  high  specific  thrust  engines  can  make  a  particularly  intense  type  of  noise 
that  has  come  to  be  called  crackle.  Its  rough,  sporadic  character  is  sometimes  likened  to  super-loud  paper 
tearing  or  the  sharp  staccato  of  water  drops  added  to  hot  oil.  It  is  unmistakably  different  than  lower  speed 
jet  noise.  Its  frequency  content,  extreme  intensity,  and  sporadic  character  make  it  particularly  annoying  and 
potentially  injurious  to  insufficiently  protected  hearing.  Observations  suggest  that  it  consists  of  series  of 
intense  shock-like  sound  waves,  which  arrive  intermittently  in  groups.  Between  these  shock  groups  a  less 
intense  but  still  loud  standard  jet  noise  is  heard.  The  steepened  shock-like  character  is  suggestive  of 
nonlinear  development,  but  to  what  degree  they  attain  this  character  at  generation  versus  nonlinear 
steepening  during  propagation  is  unclear.  There  is  evidence  that  both  are  key  factors  in  crackle.  It  is  also 
unclear  whether  the  intermittency  of  crackle  is  due  to  the  changing  character  of  the  turbulence  or  if  it  is  a 
stochastic  propagation 
effect. 

We  have  complete  a  set  of  detailed  direct  numerical  simulations  of  free  shear  flows  to  investigate 
turbulence  as  a  source  of  jet  crackle.  This  report  includes  (1)  detailed  documentation  of  these  simulations, 
(2)  their  verification  and 


validation,  (3)  and  several  post-processing  investigation  the  underlying  mechanics  of  jet  crackle,  including 
detailed  quantification  of  the  near  acoustic  field,  analysis  of  statistical  metrics  of  jet  crackle,  analysis  of 
sound-field  nonlinear  effects,  and  the  relation  of  crackle  to  the  linear  instability  modes 
supported  by  free  shear  flows.  This  is  the  first-ever  such  detailed  investigation  of  these  mechanisms.  It  has 
revealed  mechanisms  that  can  potential  be  harnessed  to  suppress  this  type  of  noise. 
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